---
title: "External Validity: Framework, Design and Analysis"
subtitle: "Main Manuscript Analyses"
author: "Naoki Egami and Erin Hartman"
date: \today
output:
  html_document:
    df_print: paged
---


```{r setup, include=FALSE}
rm(list = ls())
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(message = FALSE)
```

```{r load_libraries, include = FALSE}
## Load the evalid functions from Github.  At time of acceptance, the repository was at commit
## `f4d81f0e8dfc58f6fca8e9e3ac272bb30470cc6d` which is loaded below.

if(!require(evalid)) {
  library(devtools)
  install_github("naoki-egami/evalid@f4d81f0e8dfc58f6fca8e9e3ac272bb30470cc6d", dependencies = TRUE)
}
library(evalid)


library(grf)
library(foreign)
library(lmtest)
library(sandwich)
library(survey)
library(MASS)
library(parallel)
library(estimatr)
library(bartCause)
library(tidyverse)
library(knitr)
library(kableExtra)
library(estimatr)
```

## System Information
```{r}
print(Sys.time())

start_time <- Sys.time()

print(sessionInfo())
```

```{r}
# Generate folders for results if they do not exist
if(!file.exists("./generated/main_results")) dir.create("./generated/main_results")
if(!file.exists("./generated/Figures")) dir.create("./generated/Figures")
if(!file.exists("./generated/Tables")) dir.create("./generated/Tables")
```

```{r}
set.seed(29384)
```

## Broockman and Kalla (2018) -- Figure 7 (Main Manuscript)

``` {r}
## Load the data.

load("./generated/data/recoded_broockman_kalla_data.RData")

## Load previous results if they exist

if(file.exists("./generated/main_results/results_bk_external.RData")) {
  load("./generated/main_results/results_bk_external.RData")
}

## Set number of bootstrap replicates

nboot = 1000
```

### Overlap covariates

Estimation of individual level treatment effects (via `bartCause` and fully saturated OLS) and IPW weights (via `calibration`) is done using the overlap set of covariates across the two datasets.  SATE regression and weighted regression is done controlling for the pre-specified controls (in the original authors' PAP, and used in their final analysis).

```{r define_covariates}
## Covariates used for estimation

## Projection models: uses raw age
covariates_proj <- c('vf_age', 'vf_female', 'vf_black', 'vf_white',
                        names(sample)[str_detect(names(sample), "ideology_t0_factor|religious_t0_factor|pid_t0_factor")][-c(1:3)])

## Weighting models: uses age buckets
covariates <- c('vf_female', 'vf_black', 'vf_white',
                        names(sample)[str_detect(names(sample), "ideology_t0_factor|religious_t0_factor|pid_t0_factor|vf_age_bucket")][-c(1:4)])
  
# Pre-specified regression controls used in original analysis  
control_vars <- c('miami_trans_law_t0', 'miami_trans_law2_t0', 'therm_trans_t0', 
  'gender_norms_sexchange_t0', 'gender_norms_moral_t0', 'gender_norms_abnormal_t0',
  'ssm_t0', 'therm_obama_t0', 'therm_gay_t0','vf_democrat', 'ideology_t0', 
  'religious_t0', 'exposure_gay_t0', 'exposure_trans_t0', 'pid_t0', 'sdo_scale',
  'gender_norm_daugher_t0', 'gender_norm_looks_t0', 
  'gender_norm_rights_t0', 'therm_afams_t0', 'vf_female', 'vf_hispanic',
  'vf_black', 'vf_age', 'survey_language_es', 'cluster_level_t0_scale_mean')
```

### Trans. Tolerance Index Analysis

This section re-analyzes the initial findings from the paper, using the different PATE estimation strategies.  The estimators are:

* SATE (Regression) is the bivariate regression with pre-specified controls and HC2 errors, accounting for household-level clusters.
* PATE (IPW) is the weighting based estimator estimated using calibration weights.
* PATE (wLS) is the IPW weighted least squares, controlling for pre-specified controls, using calibration weights.
* PATE (OLS-Outcome) is OLS outcome based estimator (estimated with overlap covariates), estimated separately for Y1 and Y0, projected on to the population.  The doubly robust AIPW estimator, using calibration weights, is also included.
* PATE (BART-Outcome) is BART outcome based estimator (estimated with overlap covariates) projected on to the population.  The doubly robust AIPW estimator, using calibration weights, is also included.

**Note that the point estimates aren't identical to those presented in the original authors' paper for the DiM because they are the ITT estimates, rather than the CACE estimates.**  There are `r nboot` bootstrap replicates used to calculate the standard errors.

```{r trans_tol_index_analysis, results = 'hide', message = FALSE}
## Only re-runs analysis if no results are loaded.

## run analysis for each time period
if(!exists("plot_data_tti")) {
  plot_data_tti <- NULL
  outcome_vars <- c("trans.tolerance.dv.t1" , "trans.tolerance.dv.t2",
                    "trans.tolerance.dv.t3", "trans.tolerance.dv.t4")
  for(outcome in outcome_vars) {
    
    # limit to rows without missing outcome
    exp_data <- sample %>% filter(!is.na(sample[, outcome]))
   
    # create formulas for analysis with outcome
    formula_outcome <- as.formula(paste0(outcome, " ~ ", paste0(covariates, collapse = " + ")))
    formula_outcome_proj <- as.formula(paste0(outcome, " ~ ", paste0(covariates_proj, collapse = " + ")))
    formula_outcome_wls <- as.formula(paste0(outcome, " ~ ", paste0(control_vars, collapse = " + ")))
    formula_weights <- as.formula(paste0("S ~ ", paste0(covariates, collapse = " + ")))
    
    ## Sate Estimate
    sate_est <- data.frame(estimator = c("SATE"),  type = "SATE", 
                       tpate(formula_outcome = formula_outcome_wls,
                               formula_weights = formula_weights,
                               exp_data = exp_data, pop_data = exp_data,
                               treatment = "treat_ind", 
                               id_cluster = exp_data$block_ind,
                               est_type = "wls", weights_type = "calibration",
                               sims = nboot,
                               compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")] %>% data.frame())
    
    ## PATE Estimates
    
    ## Weighting: IPW
    ipw_cal <- data.frame(estimator = "ipw-cal",  type = "weighting", 
                       tpate(formula_outcome = formula_outcome,
                             formula_weights = formula_weights,
                             exp_data = exp_data, pop_data = pop,
                             treatment = "treat_ind", 
                             id_cluster = exp_data$block_ind,
                             est_type = "ipw", weights_type = "calibration",
                             sims = nboot,
                             weights_max = 10,
                             compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
                       %>% data.frame())
    
    ## Weighting: Weighted Least Squares
    wls_cal <- data.frame(estimator = "wls-cal",  type = "weighting", 
                       tpate(formula_outcome = formula_outcome_wls,
                             formula_weights = formula_weights,
                             exp_data = exp_data, pop_data = pop,
                             treatment = "treat_ind", 
                             id_cluster = exp_data$block_ind,
                             est_type = "wls", weights_type = "calibration",
                             sims = nboot,
                             weights_max = 10,
                             compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
                       %>% data.frame())
    
    ## Outcome: Regression
    out_ols <- data.frame(estimator = "OLS-proj", type = "outcome",
                        tpate(formula_outcome = formula_outcome_proj,
                              formula_weights = formula_weights,
                              exp_data = exp_data, pop_data = pop,
                              treatment = "treat_ind",
                              id_cluster = exp_data$block_ind,
                              est_type = "outcome-ols",
                              sims = nboot,
                              weights_max = 10,
                              compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
                        %>% data.frame())
    
    ## Outcome: BART
    out_bart <- data.frame(estimator = "BART-proj", type = "outcome",
                        tpate(formula_outcome = formula_outcome_proj,
                              formula_weights = formula_weights,
                              exp_data = exp_data, pop_data = pop,
                              treatment = "treat_ind",
                              id_cluster = exp_data$block_ind,
                              est_type = "outcome-bart",
                              sims = nboot,
                              weights_max = 10,
                              compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
                        %>% data.frame())
    ## Doubly Robust Estimators
    
    ## Doubly Robust: Augmented-Regression
    dr_cal_ols <- data.frame(estimator = "DR-OLS-cal", type = "doubly robust", 
                        tpate(formula_outcome = formula_outcome,
                              formula_weights = formula_weights,
                              exp_data = exp_data, pop_data = pop,
                              treatment = "treat_ind",
                              id_cluster = exp_data$block_ind,
                              weights_type = "calibration",
                              est_type = "dr-ols",
                              sims = nboot,
                              weights_max = 10,
                              compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
                        %>% data.frame())
    
    ## Doubly Robust: Augmented-BART
    dr_cal_bart <- data.frame(estimator = "DR-BART-cal", type = "doubly robust", 
                            tpate(formula_outcome = formula_outcome_proj,
                                  formula_weights = formula_weights,
                                  exp_data = exp_data, pop_data = pop,
                                  treatment = "treat_ind",
                                  id_cluster = exp_data$block_ind,
                                  weights_type = "calibration",
                                  est_type = "dr-bart",
                                  sims = nboot,
                                  weights_max = 10,
                                  compute_sate = FALSE)$tpate[c("est", "se", "ci_lower", "ci_upper")]
                            %>% data.frame())
    
    results = data.frame(time = which(outcome_vars == outcome),
                         bind_rows(
                            sate_est,
                            ipw_cal,
                            wls_cal,
                            out_ols,
                            out_bart,
                            dr_cal_ols,
                            dr_cal_bart
                          )
                        )
    
    plot_data_tti <- bind_rows(plot_data_tti, results)

  }
  
  
  # Format name columns
  plot_data_tti$time_name = factor(plot_data_tti$time, labels = c("+3 Days", "+3 Weeks", "+6 Weeks", "+3 Months"))
  
  plot_data_tti$estimator_type = factor(case_when(str_detect(plot_data_tti$type, "SATE") ~ "SATE",
                                         str_detect(plot_data_tti$type, "outcome") ~ "T-PATE:\nOutcome-based Estimator",
                                         str_detect(plot_data_tti$type, "doubly robust") ~ "T-PATE:\nDoubly Robust Estimator",
                                         TRUE ~ "T-PATE:\nWeighting-based Estimator"),
                                        levels = c("SATE", "T-PATE:\nWeighting-based Estimator",
                                                   "T-PATE:\nOutcome-based Estimator",
                                                   "T-PATE:\nDoubly Robust Estimator"))
  
  plot_data_tti$estimator_name = factor(case_when(plot_data_tti$estimator == "SATE" ~ "OLS",
                                                  plot_data_tti$estimator == "ipw-cal" ~ "IPW",
                                                  plot_data_tti$estimator == "wls-cal" ~ "wLS",
                                                  plot_data_tti$estimator == "OLS-proj" ~ "OLS",
                                                  plot_data_tti$estimator == "BART-proj" ~ "BART",
                                                  plot_data_tti$estimator == "DR-OLS-cal" ~ "AIPW\nwith OLS",
                                                  plot_data_tti$estimator == "DR-BART-cal" ~ "AIPW\nwith BART"),
                                                  levels = c("OLS", "IPW", "wLS", "BART", "AIPW\nwith OLS", "AIPW\nwith BART"))

  ## Save results
  save(plot_data_tti, file = paste0("./generated/main_results/results_bk_external.RData"))  
  
}
```

Generate Figure 7.
```{r plot_figure_7, echo = TRUE, fig.width = 10, fig.height = 10, fig.align='center'}
dw = 0.5

plot_data_tti %>%
  ggplot(aes(x = estimator_name, y = est, color = estimator_type)) + 
  scale_color_manual(values=c("black", "red", "blue", "magenta")) + 
  geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), position = position_dodge(width = dw)) +
  geom_point(aes(y = est), position = position_dodge(width = dw)) +
  theme_bw() +
  xlab("") +
  ylab("Estimated Causal Effects on Transgender Tolerance Scale") +
  theme(axis.title.y = element_text(size = 13, margin = margin(t = 0, r = 10, b = 0, l = 0)))+
  geom_hline(yintercept = 0, linetype = 2) +
  theme(legend.position="none",
        axis.text.x = element_text(angle = 0, size = 11, margin = margin(t = 10, r = 0, b = 0, l = 0), 
                                   colour = "black"),
        axis.text = element_text(size = 11),
        axis.title = element_text(size = 11),
        legend.text = element_text(size = 11),
        panel.grid.major = element_blank(),
        strip.text.x = element_text(size = 11),
        strip.text.y = element_text(size = 13)) +
  facet_grid(time_name ~ estimator_type, scales = "free_x")

ggsave("./generated/Figures/Figure_7.pdf", height = 8, width = 11)
```

```{r numerics_figure_7, include = FALSE, eval = TRUE}
## Outputs numeric code to "Fig7_numerics.tex"
rownames(plot_data_tti) <- NULL
cat(plot_data_tti %>%
  select(estimator_type, estimator_name, est, se, ci_lower, ci_upper) %>%
    mutate(estimator_type = as.character(estimator_type),
           estimator_type = str_replace(estimator_type, "\n", " "),
           estimator_name = as.character(estimator_name),
           estimator_name = str_replace(estimator_name, "\n", " "),
           ci = paste0("[", round(ci_lower, 3), ", ", round(ci_upper, 3), "]")) %>%
    select(-ci_lower, -ci_upper) %>%
  kable(caption = "Numeric Values for T-PATE Estimates for Broockman and Kalla (2016).",
        booktabs = T,
        col.names = c("Estimator Type", "Estimator", "Estimate", "SE", "95% CI"),
        digits = 3,
        format = "latex",
        align = c("l", "l", "r", "r", "c")) %>%
  kable_styling(latex_options = c("hold_position")) %>%
  pack_rows("Measurement at +3 Days", 1, 7, latex_gap_space = "1em") %>%
  pack_rows("Measurement at +3 Weeks", 8, 14, latex_gap_space = "1em") %>%
  pack_rows("Measurement at +6 Weeks", 15, 21, latex_gap_space = "1em") %>%
  pack_rows("Measurement at +3 Months", 22, 28, latex_gap_space = "1em"),
  file = "./generated/Tables/Fig7_numerics.tex")
```

## Bisgaard Analysis -- Figure 8 (Main Manuscript)

Load the data.
```{r Bisgaard_load_data, include = TRUE}
load("./generated/data/study1.RData")
load("./generated/data/study2.RData")
load("./generated/data/study3.RData")
load("./generated/data/study4.RData")
```

```{r Bisgaard_functions, include = TRUE}
# function for running lm
lm_analysis <- function(study, outcome_var, gov = TRUE) {
  study <- recode_study(study, outcome_var, gov)
  mod_fit <- lm_robust(outcome ~ treatment, data = study)
  return(mod_fit)
}

# recode study outcomes
recode_study <- function(study, outcome_var, is_gov_supporter) {
  if(!("resp1" %in% names(study))) {
    study$resp1 <- rep(NA, nrow(study))
  }
  if(!("c1raw" %in% names(study))) {
    study$c1raw <- rep(NA, nrow(study))
  }
  if(!("obamaspend" %in% names(study))) {
    study$obamaspend <- rep(NA, nrow(study))
    study$obamaleader <- rep(NA, nrow(study))
    study$congress <- rep(NA, nrow(study))
    study$fed <- rep(NA, nrow(study))
    study$wallst <- rep(NA, nrow(study))
    study$export <- rep(NA, nrow(study))
  }
  
  study <- study %>%
    filter(!is.na(gov)) %>%
    filter(gov == as.numeric(is_gov_supporter)) %>%
    filter(treatment %in% c("pos", "neg")) %>%
    mutate(resp1 = case_when(is.na(resp1) ~ 0.5,
                             TRUE ~ as.numeric(resp1)),
           c1 = case_when(c1raw == 0 ~ 1,
                          c1raw == 3 ~ 2,
                          c1raw == 15 ~ 2,
                          c1raw == -1 ~ 1,
                          TRUE ~ as.numeric(c1raw)) - 1,
           obamaspend = obamaspend - 4,
           obamaleader = obamaleader - 4,
           congress = -1 * (congress - 4),
           fed = -1 * (fed - 4),
           wallst = -1 * (wallst - 4),
           export = -1 * (export - 4)) %>%
    mutate(treatment = as.numeric(treatment == "neg"),
           outcome = get(outcome_var))
  return(study)
}
```

Run analysis for opposition supporters.
```{r Bisgaard_opposition_analysis_opp, include = TRUE}
## For Opposition Supporters
opp_results <- list(
    # (1) Close, Study 1 in the US 
    lm_analysis(study1, "resp1", 0),
    # (2) Close, Study 2 in Denmark 
    lm_analysis(study2, "resp1", 0),
    # (3) Close, Study 3 in Denmark 
    lm_analysis(study3, "resp1", 0),
    # (4) Open, Study 1 in the US 
    lm_analysis(study1, "c1", 0),
    # (5) Open, Study 3 in Denmark 
    lm_analysis(study3, "c1", 0),
    # (6) Open, Study 4 in Denmark
    lm_analysis(study4, "c1", 0),
    # (7) Argument/Pro -- obamaspend, Study 1 in the US
    lm_analysis(study1, "obamaspend", 0),
    # (8) Argument/Pro -- obamaleader, Study 1 in the US
    lm_analysis(study1, "obamaleader", 0),
    # (9) Argument/Con -- congress, Study 1 in the US
    lm_analysis(study1, "congress", 0),
    # (10) Argument/Con -- fed, Study 1 in the US
    lm_analysis(study1, "fed", 0),
    # (11) Argument/Con -- wallst, Study 1 in the US
    lm_analysis(study1, "wallst", 0),
    # (12) Argument/Con -- export, Study 1 in the US
    lm_analysis(study1, "export", 0)
)

opp_results <- lapply(opp_results, function(x) (tidy(x) %>% filter(term == "treatment"))[, c("estimate", "std.error")]) %>% bind_rows()

opp_pct <- data.frame(type = "(H2) Opposition Supporters", pct(1 - pnorm(opp_results$estimate/opp_results$std.error)))
```

Run analysis for government supporters.
```{r Bisgaard_opposition_analysis_gov, include = TRUE}
## For Government Supporters
gov_results <- list(
    # (1) Close, Study 1 in the US 
    lm_analysis(study1, "resp1", 1),
    # (2) Close, Study 2 in Denmark 
    lm_analysis(study2, "resp1", 1),
    # (3) Close, Study 3 in Denmark 
    lm_analysis(study3, "resp1", 1),
    # (4) Open, Study 1 in the US 
    lm_analysis(study1, "c1", 1),
    # (5) Open, Study 3 in Denmark 
    lm_analysis(study3, "c1", 1),
    # (6) Open, Study 4 in Denmark
    lm_analysis(study4, "c1", 1),
    # (7) Argument/Pro -- obamaspend, Study 1 in the US
    lm_analysis(study1, "obamaspend", 1),
    # (8) Argument/Pro -- obamaleader, Study 1 in the US
    lm_analysis(study1, "obamaleader", 1),
    # (9) Argument/Con -- congress, Study 1 in the US
    lm_analysis(study1, "congress", 1),
    # (10) Argument/Con -- fed, Study 1 in the US
    lm_analysis(study1, "fed", 1),
    # (11) Argument/Con -- wallst, Study 1 in the US
    lm_analysis(study1, "wallst", 1),
    # (12) Argument/Con -- export, Study 1 in the US
    lm_analysis(study1, "export", 1)
)

gov_results <- lapply(gov_results, function(x) (tidy(x) %>% filter(term == "treatment"))[, c("estimate", "std.error")]) %>% bind_rows()

gov_pct <- data.frame(type = "(H1) Government Supporters", pct(pnorm(gov_results$estimate/gov_results$std.error)))
```

Combine and code the type of outcome question and context.
```{r bisgaard_color_results, include = TRUE}
res <- bind_rows(opp_pct, gov_pct) %>%
  mutate(Yval = case_when(
    h_num %in% c(1, 2, 3) ~ "Open-Ended",
    h_num %in% c(4, 5, 6) ~ "Close-Ended",
    TRUE ~ "Argument Rating"),
    Cval = case_when(
      h_num %in% c(1, 4, 7, 8, 9, 10, 11, 12) ~ "United States",
      TRUE ~ "Denmark"
    ))

# Rename
res$type = as.character(res$type)
res$type[res$type == "(H1) Government Supporters"] <- "(H1) Incumbent Supporters"
res <- res[c(which(res$type == "(H1) Incumbent Supporters"), which(res$type != "(H1) Incumbent Supporters")), ]

save(res, file = paste0("./generated/main_results/bisgaard_res.RData")) 
```

Generate Figure 8.
```{r include = TRUE}
res_plot <- ggplot(res, aes(x = threshold, y = p_value, color = Cval, shape = Yval)) + geom_point() + 
  ylim(c(0, 0.3)) + geom_hline(yintercept = 0.05, linetype = 2) + 
  theme_bw() +
  xlab("Threshold") + ylab("Partial Conjunction p-value") +
  theme(axis.title.y = element_text(size = 10, margin = margin(t = 0, r = 10, b = 0, l = 0)),
        axis.title.x = element_text(size = 10, margin = margin(t = 10, r = 0, b = 0, l = 0)),
        legend.title=element_text(size=9)) +
  facet_wrap(~type) +
  # geom_text(x=2, y=0.28, label="8/12") + 
  scale_color_manual(guide = guide_legend(title = "Context"), values=c("red", "blue")) +
  scale_shape(guide = guide_legend(title = "Outcome")) +
  scale_x_continuous(breaks = c(2, 4, 6, 8, 10, 12)) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank())

res_plot
ggsave("./generated/Figures/Figure_8.pdf", width = 6, height = 2.5)
```

```{r numerics_figure_8, include = FALSE, eval = TRUE}
## Outputs numeric code to "numerics_main.tex"

row.names(res) <- NULL
cat( "\n\n\n", 
res %>%
       select(threshold, Yval, Cval, p_value) %>%
     kable(caption = "Numeric Values for Sign-Generalization Test for Bisgaard (2019).",
        booktabs = T,
        col.names = c("Threshold", "Outcome Variation", "Context Variation", "P-value"),
        digits = 3,
        format = "latex") %>%
  #kable_styling(latex_options = c("striped")) %>%
  pack_rows("(H1) Incumbent Supporters", 1, 12, latex_gap_space = "1em") %>%
  pack_rows("(H2) Opposition Supporters", 13, 24, latex_gap_space = "1em"),
     file = "./generated/Tables/Fig8_numerics.tex")
```

Total runtime is `r round(difftime(Sys.time(), start_time, units = "mins"), 2)` minutes.
